Iranian Journal of Numerical Analysis and Optimization 
Vol. 9, No. 2, (2019), pp 123-139 Poxso 


DOI:10.22067 /ijnao.v9i2.75991 yu 


Research Article 


An extension of the quasi-Newton 
method for minimizing locally 
Lipschitz functions 


Z. Akbari * 


Abstract 


We present a method to minimize locally Lipschitz functions. At first, a 
local quadratic model is developed to approximate a locally Lipschitz func- 
tion. This model is constructed by using the e-subdifferential. We minimize 
this local model and compute a search direction. It is shown that this direc- 
tion is descent. We generalize the Wolfe conditions for finding an adequate 
step length along this direction. Next, the method is equipped with a quasi- 
Newton approach to update the local model and its globally convergence 
is proposed. Finally, the proposed algorithm is implemented in MATLAB 
environment on some standard nonsmooth optimization test problems and 
compared with some algorithms in the literature. 
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1 Introduction 


In this paper, we consider the following unconstrained nonsmooth optimiza- 
tion problem: 

i 1 

min f(x) (1) 

where f : R” — R is a locally Lipschitz continuous function. There ex- 


ist several methods for solving this problem, for example, the subgradient 
method [25], the bundle methods [7,10,19,24], algorithms based on smoothing 
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techniques [22], the derivative free methods [3], the bundle and quasi-Newton 
combining methods [13,15,17], the gradient sampling method [5,12], the lim- 
ited memory gradient bundle method [11], and the trust region method [2,23]. 


In the most smooth optimization problems, the objective function is ap- 
proximated by a quadratic model as follows: 


, d’ B,d 

Qlen +d) = flee) + VF(oe)7d + (2) 
where V f(a.) and By, are the gradient and Hessian of f at x,; see [21]. In 
fact, the objective function is locally approximated at x,. In quasi-Newton 
methods, the quadratic model (2) is minimized and its solution is a descent 
direction at x,. Afterward a line search method is applied along this direc- 
tion. This process is repeated until the optimal solution is achieved. Also, 
in trust region methods, the model (2) is minimized on a given region, called 
trust region, and a search direction is obtained. If this direction does not 
decrease the objective function adequately, then the trust region is updated. 


In the nonsmooth optimization same as the smooth one, the objective 
function is approximated. Han et al. [9], presented a quadratic model to 
approximate the locally Lipschitz function f as follows: 


d’ Byd 
9 ? 


Q(ae +d) = f(we) + b(ax, d) + (3) 
where ¢: R” x R" — R is a given iteration function and By is a given n x n 
symmetric matrix. The iteration function ¢ must have some properties to 
guarantee the convergence of produced minimization algorithm. 


The presented iteration function in [9] is not practical. So, in this paper, 
we propose a practical iteration function and solve the following problem: 


d’ Byd 


min Qe(xn + d) = f(x) + (xe, d) + —J—- (4) 


Suppose that d, is the solution of (4). We show that dy is a descent direction 
and develop a line search method to reduce objective function along this 
direction. 


In this paper, we propose a computable quadratic model to approximate 
the locally Lipschitz function and develop a minimization algorithm based 
on this model. In Section 2, some preliminary concepts of the nonsmooth 
analysis are reviewed. In Section 3, an iteration function is introduced for 
the locally Lipschitz objective function. After that, a descent algorithm is 
proposed. In Section 4, the global convergence of the presented algorithm 
is proved. The line search algorithm is described in Section 5. Numerical 
results are given in Section 6. 


An extension of the quasi-Newton method for minimizing ... 125 


2 Preliminaries 


In this section, we state the basic concepts and definitions of nonsmooth 
analysis from [4]. The Clarke generalized directional derivative of the locally 
Lipschitz function f at the point x in the direction d is defined by 


f° (a, d) := limsup fly + td) — f(z) 
— ya,tl0 t . 


The Clarke generalized gradient of f at x is the set Of (x) defined by 
Of (z) = {€ € R"| f° (a; d) > €7d for all d€ R"}. 


Each vector € € Of (x) is called a subgradient of f at 2. For € > 0, The Clarke 
generalized e-directional derivative of f at x in the direction d is defined by 


Fe (a, d) = limsup f(y + td) — f(z) + 
‘ ya, tlo t ’ 


and the Goldstein ¢-subdifferential of f at the point x is the set 


d.f(x) == cl conv{Of(y), |]a—ylle <4, 


where conv and cl are used to denote the convex hull and the closure, respec- 
tively. Each element € € 0. f(a) is called an e-subgradient of the function f 
at x; see [4]. It can be seen that f2(x,d) = supecg, f(x) é?d. Let L be the 
Lipschitz constant of f in a neighborhood of x. Then, we have 


lvll2< £, for all v € Of (a). 


Throughout the paper, the used norm is the Euclidean norm. If f is dif- 
ferentiable at x, then Vf(x) € Of(x). Furthermore, if f is continuously 
differentiable at x, then 


Of (x) = {VF(x)}- 


If 0 € 0. f(x), then x is called as an e-stationary point. 


3 Generalized quasi-Newton method 


In this section, we develop a computable iteration function and suppose that 
By, is positive definite and has bounded norm, m < || Bx|| < M. Since B,° 
is positive definite, then we consider its Cholesky decomposition as follows: 


Boat b, 
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where L is an upper triangular matrix with real and positive diagonal entries, 
and L? denotes the transpose of L. To define the iteration function, we 
consider the following problem: 


Vp =arg min ||Lé|l’, 5 
p=arg min || (5 


where ¢ is a small enough positive scalar, and define ¢(x;,d) = vjd. There- 
fore the problem (4) is rewritten as follows: 


, d? Bid 
min qi(d) = upd + aa: (6) 


By this modification, —By ug is the unique solution of (6). Let dy, = 
By up 
(By eelle 
rection and the step length a, is returned. The next iteration is computed 


as follows: 


be a search direction. The line search is applied along this di- 


Le+1 = Te + aed. 


At the point 7,41, either B;, is updated by one of the quasi-Newton methods 
or is set the identity matrix for all k. Now, we are ready to express the 
generalized quasi-Newton algorithm for solving the problem (1): 


Algorithm 1 Generalized quasi-Newton algorithm 


Step 0 (Initialization): Let By = Inyn, € >0, ec: € (0,1), 71 ER", k= 
1. 


Step 1 (Creating an iteration function): Set Las an upper triangular 
matrix of Cholesky decomposition for the matrix B,'. Solve (5) at the 
point x, and denote its solution by vz. 


Step 2 (Stopping condition): If ||v;,||2 = 0, then stop, else set d, = 
_ By tur 
Bz vell 


and compute gx = qx(dx). 


Step 3 (Line search method): Apply the line search algorithm and find 
the step length a, € [e, 1] such that 


Qy, = Max {a | f(a, t+ady) — f(x) < cau; dj} . 
Set Tet41 = e+ Apdg. 


Step 4 (Updating): If it is necessary, then update By, set k = k+1 and 
go to step 1. 
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As mentioned before, we can set By = I for all k in Step 4, which means 
that the matrix Be is a fixed and identity matrix. Therefore, Algorithm 
1 is converted to the steepest descent method for minimizing the locally 
Lipschitz function [18]. If Bj, is updated by the BFGS formula (see [21]), then 
Algorithm 1 is a generalization of the quasi-Newton method for minimizing 
the locally Lipschitz continuous function. In Section 5, we describe a line 
search algorithm how to find a step length along the search direction such 
that the Wolfe conditions are satisfied. 

The following lemma shows that the optimal value of (6) is nonpositive. 
Also it ensures that if the optimal value of (6) is nonzero, then the corre- 
sponding solution provides a descent direction for the function f at the point 
xr;,. Else x, is the optimal solution. 


Lemma 1. Let f : R" > R be a locally Lipschitz function and let c, € (0,1). 
Suppose that vy, is the solution of problem (5) at x, such that vy #0. If 


d, = By ‘ve then 
kB, vella’ 


(a) vp dx <0 and q, <0. 


(b) if upd, > 0, then there exists a scalar @ > 0 such that for all a € [0, a], 
we have 
f (ae + adg) — f (ag) < cravz dg. 


Proof. Part (a) is clear. To prove part (b), at the contrary, we suppose there 
exists a positive sequence {a;} such that a; — 0 and we have 


f(z a aidy) — f (ar) ea c1ayup dp. 


Both sides of the above inequality are divided by a;, then, with passing to 
the lim sup,_,,,, we deduce 


f(x + idk) — (ak) 


fe (xk, dk) = f° (we, de) = lim sup > cup de. 
ioo ay 
_ __ By lve Tq, — — eB ve — __[Lvell3 
Since d, = [B-Tvele? then uv; dp = [Botella ~~ 1Betonla” Thus we have 
. || Lv. ||? 
je See) 2G (7) 
|By *velle 
On the other hand, we have 
i T Boi 
fe(an,dx) = max €"d, = max ae 
£0. f (ex) €€0- f(r) ||B, velle 
—1 —1 
in €7L7 Ly, = —~—— in (LE)" Lug. 


~ TB-u. lle cearr ~ 7 B-t = 
| By “Ux|l2 Oe f(ex) | Bol uglle €€8. (ex) 
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According to (5), we have ||Lv,||? < ||LE||?, for all € € 0,f(x), and since the 
set O.f(a,) is a convex compact set, referring to [20, Lemma 5.2.6.], reader 
can easily find out 


° || Len ll3 
Lp, dy) = —-———. 8 
tel k k) [Bz valle ( ) 
By (7) and (8), we have 
Lvxl2 ; Lv,||2 
_ | vel = f (xz, de > —cy | ella : 
|B elle | Ba velle 


Since c; € (0,1), the above equation is false. Thus, part (b) is proved. 


The following lemma shows that the step length a, = € is satisfied the 
Arimijo condition along d, for all k. 


Lemma 2. Let f : RR” > R be a locally Lipschitz continuous function and 
let c, € (0,1). Then the step length € satisfies the Armijo condition along dy, 
at Xp. 


Proof. Since f is a locally Lipschitz continuous function, then according to 
the mean-value theorem in [20], there exists 7 € (0,1) such that z = x, +ned, 
and 

f (ay + ed) — f(xy) € OF (z)7 (Edy). 


So, there exists u € Of(z) such that 
eu? Br up (Lu)? (Luz) 


f (ax + edy) — f (wx) = eu™ dy = = = = 
Cte k) (xx) k |Bztvelle Bz velle 


Since vz € O-f (xz) is the solution of problem (5), and 0.f (x) is a convex 
compact set, then by [20, Lemma 5.2.6], we have 


(Lu)? (Lux) > ||Loxl3- (10) 
So, by (9) and (10), we have 


|| Lvell5 cup By Up T 
f(a, + edz) — flap) < = = = = ev, dy. 
|B ‘velle By 'velle ‘ 


Since upd, <0 and c; € (0,1), then 


f (ap + edy) — f (ap) < crevp dr. 


Thus the proof is completed. 
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4 Optimality condition and global convergence 


In this section, we provide the necessary optimality condition and prove the 
global convergence of Algorithm 1. If ||v,||2 = 0 at some iteration k, then 
0 € O.f(a,). This shows that x, is an e-stationary point. Now, we show 
that if Algorithm 1 dose not terminate after finite many iterations, then each 
accumulation point of the sequence {x,} is an estationary point of f. To 
prove this, we first prove the following lemma. 


Lemma 3. Suppose that the level set L:= {x : f(x) < f(ai)} ts bounded, and 
there exist positive constants M and m such that ml||d||3 < d7 Byd < M|\d\\3, 
for all k and d € R”. If Algorithm 1 does not terminate after finite many 
iterations, then 

lim ||vx|| = 0. 

k—-00 


Proof. Suppose that Algorithm 1 does not terminate after finite many itera- 
tions. Since ||v,|| 4 0 for each k, then by Lemma 1, we have 


f(x + andy) — f(ar) < cragug dr. (11) 


Since {f(x,)} is a strictly decreasing and bounded below sequence, then it 
converges and we have 


{te Onde) — Fee) > 0. 


So by (11), we have 
0<c, lim pvp dp. 
k-> oo 


Tp-l 
v;, By Vr 


By Lemma 1, we know vf d, = —=+—* < 0, so 
|B," velle 
. Z _ agup By uR 
O<c, lim a,zv,z dy, = —c, lim a < 
k- oo k—- 00 |B; Up|l2 


On the other hand, by Lemma 2, a, > € Thus 


Tp-l 
B 
lit pct me =, (12) 
k-+o0 || Bi, vel 


Let 0 < \¥ < MK < -- < \* be the eigenvalues of B,*. Since B,;' is a 
positive definite matrix, then we have \¥ > m and \* < M. Therefore 


m ei vt Br up Nui 
qyllvell2 S Tr lleslle < [Bonde = ; 
n 


M 
lt —= : 13 
Boul lea SF llvelle (13) 


Then, by (12) and (13), we conclude ||v,||2 — 0. 
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In the following theorem, we show that each accumulation point of the 
sequence {x;,,} is an e-stationary point. 


Theorem 1. Suppose that the level set L:= {x : f(x) < f(x1)} ts bounded 
and that f : R” > R is a locally Lipschitz function on L. If x* is an 
accumulation point of the sequence {xz}, then x* is an €-stationary point of 


f. 


Proof. Since {x,} C L, then the sequence {z,} is bounded and it has a 
convergent subsequence. Let {x,,} be a subsequence of the sequence {2,} 
convergent to x*. We have ux, € Of (ax,) and by Lemma 3, vx, > 0. On 
the other hand, 0.f(-) is upper semicontinuous. Thus 0 € 0.f(2*) and this 
shows that x* is an e-stationary point of f. 


5 Line search algorithm 


In this section, we present a line search algorithm to update the matrix By 
by the BFGS method. To do this, we need two subgradients € € Of(a,) and 
€e41 © Of (a@p41), such that the curvature condition is satisfied, that is, 


Eky ide 2 coh di, 


where cg € (c1,1). To find £441, we use the line search strategy same as 
[21,26]. The step length a must satisfy the Wolfe conditions at x, along dx. 
The Wolfe conditions are the Armijo and curvature conditions, as 


f(a, + ad,) — f(ap) < cafe (xp, dy), (Armijo condition) 
Exa1 € Of (xe + adg), 8.t. &h41dx > cof? (xe, dx). (curvature condition) 


In this paper, we present an algorithm to find a step length, which satisfies the 
Armijo condition. Now we present the following algorithm. This algorithm 
finds an interval including a step length satisfies the Wolfe conditions at the 
point x, along d,. Algorithm 2 is started with ag = €, because the step 
length € satisfies the Armijo condition. The following proposition shows that 
Algorithm 2 terminates after finite many iterations. 


Proposition 1. If f2(ax~,dz) < 0 and f is bounded below at x, along dx, 
then Algorithm 2 terminates after finite many iterations. 


Proof. Since ¢(a@) = f(a, + ad,) is bounded below and the line l(a) = 
f(rr) + craf? (xp, dx) is unbounded below, then there exists @ such that 
(a) > Ua) for alla > @. Thus Algorithm 2 terminates when a; > @. 


Either Algorithm 2 finds a step length that it satisfies the Wolfe conditions 
or it returns an interval. We prove that there exists a step length in this 
interval, which satisfies the Wolfe conditions. 
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Algorithm 2 Line search algorithm 


set ay =l,ag=6 
repeat 
if f(v, +aidx) — f(ar) > aif? (ap, dy) then 
a*=wolfe(ag, a1) 
return a* 
else if Jé € Of (xp +aidy) s.t. Edy > cof? (ap, dy) then 
aX =a, 
return a* 
end if 
ag = 1; 
ay, = 2a, 
until 


Proposition 2. In Algorithm 2, suppose that f(a, + aidx) — f(ap) > 
cai fe(ap, dp). Then the interval (ao,a1) contains step lengths satisfying 
the Wolfe conditions. 


Proof. Let (a) = f(a, + ady) — f(a) — coafe (ap, dy). It is obvious that 
wW(ao) < w(ay). Since the curvature condition is not satisfied at ao, we have 


ET dy < cof? f (ap, dy), VE € OF (ay + Aod,). 


Thus 0 ¢ 0¢(a0). Therefore ap is not the local minimizer of w on [ao, a1]. 
So the minimum point of ~ must be in (ao, a1). Suppose that a* € (ao, a1) 
is the minimizer of ~. Thus 0 € 0¢(a*). Since 


Ob(a*) C Of (ze + a* dp) de + cof? (te, dk), 


there exists € € Of (a, +a*d,) such that €7d, — cof? (xx, dp) = 0. On the 
other hand, w(ao) > w(a*) and apg satisfies the Armijo condition. Thus 


a* satisfies the Armijo condition. Therefore, a* is satisfied in the Wolfe 
conditions. 


Now we present an algorithm to find a step length, which satisfies the 
Wolfe conditions on [ap, ai]. The following proposition indicates the conver- 
gence of Algorithm 3. 


Proposition 3. If Algorithm 8 does not terminate after finite many itera- 
tions, then it converges to a* such that 0 € Of (xy +a*dz) dy — cof? (xp, dy). 


Proof. Let ~ be the function defined in Proposition 2. Since 7(a;) < w(b;) 
and 0 ¢ Of (a, + aidy) dy — cof? (xp, dy), then 0 ¢ OW(a;). Thus, 7) takes its 
minimum on (a;,b;). Suppose r; is the minimum of w on (a;, b;). So, we have 
0 € Ouw(r;). On the other hand, {a;} and {b;} are monotone and bounded 
sequences. Therefore, these sequences are convergent. Also, we have 
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Algorithm 3 Wolfe Algorithm 


t=1 

ay = Ay 

b; = a0 
— aithi 

t; = 4 5 

repeat 


if f (rr +r tidx) — f (xr) = citi fe (ar, dr) then 
1 = tj and bj44 = b; 
else if 4é € Of (x, + tidy) such that €7d, > cof? (xp, dy) then 


return {¢; 
else 
bi44 = ti and QAji+1 = AY. 
end if 
t, = withi 
a 2 
t=i+1 
until 


lim ay — b; = 0, 
ICO 
thus 
lim a; = lim }; = lim r; = a*. 
1-00 41-00 41 OO 
Since 02)(-) is upper semicontinuous and 0 € Ow(r;), then 0 € Ov(a*). This 
shows that 0 € Of (a, + a*dg)? dy — cof? (arp, dy). 


Now, we go back to problem (5). Solving this problem is impractical, 
because the structural of the set 0,f(x) is unclear. So the set 0.f(x) is 
approximated and problem (5) is approximately solved. In [1], the generalized 
h-increasing point algorithm is proposed for computing an approximation 
solution of problem (5). Here, we summarize the generalized h-increasing 
point algorithm. Let W, be a finite subset of 0. f(x,). The convex hull of W, 
is considered as an approximation of 0. f(x,). We solve the following problem 
instead of (5): 


min v? Baw, 


veconvW; 
and consider w, as its solution. Let g, = — Twa: If the inequality 
f(& + €gn) — fee) < —erellwell, (14) 


is satisfied, then convW;, is an acceptable approximation of 0. f(x,). In this 
case, we consider g, and —||w,|| as an approximation of d, and f° (xx, dx). 
Else a new element from 0.f(x,), such as vj, must be added into W; such 
that v, ¢ convW;,. To find such an element, we apply the generalized h- 
increasing point algorithm in [1]. In [1], it is proved that this procedure is 
terminated after finite many iterations. Either the direction g, is returned 
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such that (14) is satisfied or we have 
Il well < 4, (15) 


where 6 is a predefined threshold. If (14) is satisfied, then the Armijo con- 
dition is satisfied along direction g, at x, with step length «. Also when 
(15) holds, we can assume that x; is an e-stationary point. Finding such a 
direction is given in Algorithm 3.1 in [1]. 


6 Numerical experiments 


In the numerical experiments, as computing d; by (5) is impractical, we use 
Algorithm 3.1 in [1] and approximate it with g,. In this case, f2 (ax, dx) 
is approximated with —||w;||. We set parameters as follows: « = 107°, 
5 = 10-8, «| = 107+, and cp = 0.9. In Algorithms 2 and 3, the curvature 
condition is checked just by one subgradient from 0, f(x, + td,), say €, that 
is, the following condition is checked: 


ET dy > cof? (an, de). 


In this section, we used some test problems in [8,14] and report the com- 
putational and numerical results of Algorithm 1 denoted by ”NQSN”. Also 
Algorithm 1 is compared with some algorithms existing in the literatures. 
Some of the selected algorithms are implemented in MATLAB environment 
and other ones in Fortran. The measurement of algorithm efficiency is the 
number of function evaluations. To have a better comparison of the imple- 
mented algorithms, the performance profile of Dolan and More in [6] is used. 


Two classes of test functions are applied. The first class is taken from [8] 
and the second one is the TEST 29 in [14]. The test problems are introduced 
in Table 1. We compare the smooth BFGS method [21], the variable metric 
bundle method (PVAR) [16,17], a method presented in [18] (MY), the limited- 
memory BFGS method (LBFGS) [21], the limited memory bundle method 
(LMBM) [8], and the gradient sampling method (GS) with the presented 
method. 


Each algorithm is terminated after 10000 iterations. In the performance 
profile, an algorithm solves a problem successfully when 
[Fri i 7" | < thr, 
Maly 
where thr € (0,1) is a predefined threshold and fmin and f* are the optimal 
value and achieved optimal value, respectively. In this paper, we set thr equal 
to 1074. 
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Table 1: Test problems and their optimal value for n = 1000. 


First class of problems Second class of problems 
No. problem optimal value No. problem optimal value 
Hb MAXQ 0 abe problem 2 from TEST29 0 
2 MAXHILB 0 12 problem 5 from TEST29 0 
3 LQ -1.412799e+003 13 problem 6 from TEST29 0 
4 CB3I 1998 14 problem 11 from TEST29 = 1.203128e+004 
5 CB3Il 1998 15 problem 13 from TEST29 —5.661313e+002 
6 NACTFACES 0 16 problem 17from TEST29 0 
7 Brown 2 0 17 problem 19 from TEST29 0 
8 Mifflin 2 -7.065034e+002 18 problem 20 from TEST29 0 
9 Crescent I 0 19 problem 22 from TEST29 0 
10 Crescent II 0 20 problem 24 from TEST29 0 


Let n be the size of the test problem. The algorithms are tested with 3 
sizes of test problems, the small size n = 10, middle size n = 100, and large 
scale n = 1000. Since the GS algorithm is very time consuming for large 
scale problems, then we do not run this algorithm for large scale problems. 
In Figures 1 and 2, we report the results for the first and second class of test 
problems as the performance profile, respectively. 


Figure 1 shows that the first class of the test problems is simpler than the 
second one. The proposed algorithm is more efficient than other algorithms 
for middle and large scale problems and is the only algorithm can solve all 
these problems. Also, the number of function evaluations of NQSN method 
is less than other methods. In the second class of test problems, for small size 
problems, the GS method is more efficient than other ones. But the BFGS 
and NQSN have similar results. But in the middle size, the MY method can 
just solve all problems. Performance profiles show that the NQSN method 
solves problems with significantly less number of function evaluations. For 
large scale test problems, just MY and NQSN can solve 6 problems. While 
the number of function evaluations in the NQSN method is very less than 
the number of function evaluations in MY method. The results show that 
the NQSN method is more efficient than other algorithms for large scale 
problems. Also, this algorithm solves problems with less number of function 
evaluations. 
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(a). Numerical results for n = 10. 
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(b). Numerical results for n = 100. 
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(c). Numerical results for n = 1000. 


Figure 1: The performance profiles of the first class of test problems 
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(c). Numerical results for n = 1000. 
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Figure 2: The performance profiles of the second class of test problems 
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7 Conclusion 


In this paper, the generalized quasi-Newton algorithm is presented for mini- 
mizing locally Lipschitz functions. Same as the iterative optimization meth- 
ods, first we introduced a descent direction for the locally Lipschitz function 
in each iteration. Then, we computed an adequate step length along this 
direction satisfies the generalized Wolfe conditions. After that, we showed 
that the generalized quasi-Newton algorithm is convergent. The presented 
algorithm was applied to some standard test functions and compared with 
the MY method, the smooth BFGS method, the gradient sampling method, 
the limited-memory BFGS method and the limited memory bundle method. 
In future work, we will develop this algorithm for minimizing locally Lipschitz 
functions with nonlinear constraints. 
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